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Abstract 

An array of ten broadband stations was installed on the Popocatepetl volcano (Mex- 
ico) for five months between October 2002 and February 2003. 26 regional and 
teleseismic earthquakes were selected and filtered in the frequency time domain to 
extract the fundamental mode of the Rayleigh wave. The average dispersion curve 
was obtained in two steps. Firstly, phase velocities were measured in the period 
range [2 - 50] s from the phase difference between pairs of stations, using Wiener 
filtering. Secondly, the average dispersion curve was calculated by combining obser- 
vations from all events in order to reduce diffraction effects. The inversion of the 
mean phase velocity yielded a crustal model for the volcano which is consistent with 
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previous models of the Mexican Volcanic Belt. The overall crustal structure beneath 
Popocatepetl is therefore not different from the surrounding area and the velocities 
in the lower crust are confirmed to be relatively low. Lateral variations of the struc- 
ture were also investigated by dividing the network into four parts and by applying 
the same procedure to each sub-array. No well defined anomalies appeared for the 
two sub-arrays for which it was possible to measure a dispersion curve. However, 
dispersion curves associated with individual events reveal important diffraction for 
6 s to 12 s periods which could correspond to strong lateral variations at 5 to 10 
km depth. 

Key words: Volcano seismology, Popocatepetl volcano, Rayleigh waves, Crustal 
structure 



1 Introduction 



Popocatepetl is a large andesitic strato- volcano, located 60 km south-east of 
Mexico City and 40 km West of Puebla (fig. l.a). It belongs to the Trans- 
Mexican Volcanic Belt (MVB). Its large cone is the second highest summit of 
Mexico (5452 m above sea level) with an elipsoidal 600-800 m wide crater. 
The present active period began on December 21st 1994. Since 1996, an an- 
desitic to dacitic dome cyclic ly grows into the crater and bu r sts producing 



high plumes of gas and ash (jArcieniega-Ceballos et al 



2000; 



Wright et al. 



20021 ) . More than 100 000 persons could potentially be directly affected by 



an eruption and ashes could affect an area with more than 20 million people 



( iDe La Cruz-Reyna and Siebe 



1997 



Macias and Siebe 



20051 ). 
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'he overall crustal structure beneath the MVB is re 



( iCampillo et al 



1996 



Valdes et al 



1986 



Shapiro et al. 



atively well studied 



19971 ). On the con- 



trary, the crustal seismic stru cture beneath Popo c atepe tl is not well known. 



Receiver functions analysis by 



Cruz-Atienza et al. 



(120011 ). using 4 events from 



South America, indicates that a Low Velocity Zone may be present beneath 
a station located 5 km north of the crater. 

The aim of this paper is to improve the knowledge of this complex volcano 
structure, and particularly to determine if the whole crust beneath the vol- 
cano is significantly different from the rest of the MVB. The first kilome t ers o 



crust beneath several volcanoes have been studied (e.s. 



Laigle et al. 



Dawson et al. 



1999 



2000 



Benz et al 



19961 ). Typical volcanic anomalies are low ve- 
locity zones, attributed to the presence of partial melt, or high velocity zones, 
due to solidified magmatic intrusions. 

We concentrate on the S-wave structure, as S-wave velocities are very sensitive 
to temperature changes and to the presence of even small amounts of partial 
melt. The easiest way to get an overall picture of S-wave velocities is through 
surface wave analysis. However, the traditional 2-stations methods can not be 
used in this rather diffractive environnement as measurements would possi- 



bly be strongly biased due to local and regional diffraction (IWielandt 



1993 



Friederich et al. 



19951 ). An alternative approach is therefore to use array anal- 



ysis. Suc h methods have been used on volcanoes, particu larly for tremor source 



location (IMetaxian et all 120021; lAlmendros et al. 



tu re study 



111 



Chouet 



( See fo r example 



Saccorotti et al 



2002j) or for shallow struc- 



200lh . More details can be found 



(120031 ) who presents a state of the art on volcano seismology. 
The dispersion curve had to be measured over a wide frequency range (0.02-1 
Hz) to study the overall crustal structure beneath Popocatepetl. The array 



configuration which was strongly influenced by topography and logistic is- 
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sues was such that we could not use spatial Fourier transforms outside a very 
narrow frequency range. Consequently methods based on wavenumber decom- 
position were excluded. The use of time-domain methods was problematic as 
we needed a good frequency resolution. 
These considerations led us to use the procedure of 



Pedersen et al. 



((20031) to 



measure phase velocities across the array. The assumption behind this method 
is that the records are constituted by one single plane wave which propagates 
through the array. Even though this hypothesis is most probably wrong for 
most individual events, it may be corrected by averaging out unwanted waves 
(diffraction effects, non plane waves, etc) using events from different directions. 
The variability between different events will also provide an error estimate on 
the dispersion curve. To increase frequency range and azimuthal coverage we 
used both teleseismic and local events. 

After a short description of the data and the processing methods used, we 
present and discuss the main results, with a comparison of the overall crustal 
structure beneath the volcano to that of the MVB. 



2 Data 



An array of nine stations (Guralp CMG 40T) with three-components broad- 
band sensors (30-60 s cut-off period) was installed in October 2002 on the 
Popocatepetl volcano and continuously recorded four months of seismic events. 
Figure l.b shows the array geometry. The station altitudes were between 2500 
and 4300 m above sea level. The reference altitude used in this study corre- 
sponds to an average level of 3500 m a.s.l. 
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To obtain dispersion curves in a period range of 2 to 50 s, we chose to use both 
teleseimic and local events with epicentral distances between 200 and 15000 
km (see fig. 2). We selected vertical components of events with a good signal 
to noise ratio and with well developed Rayleigh waves. The usuable frequency 
range for the two types of events overlapped, however the long period part of 
the dispersion curve was mainly calculated using teleseismic events while the 
shorter periods were dominated by regional events. 

Prior to the array analysis, we deconvolved the data with the instrument re- 
sponses. The second step of this anal ysis was to enhance the signal-to-noise 



ratio through time-frequency filtering (jLevshin et al. 



19891 ). In this part of the 



analysis we firstly applied multiple filter to the data and identified the group 
velocity dispersion curve by the maximum amplitude at each frequency. We 
secondly integrated this curve to obtain the phase velocities and subtracted 
the corresponding phase 6(f) at each frequency to obtain a non dispersive 
signal. A time window was then applied onto the non- dispersive wave to sup- 
press noise and the phase 6(f) was finally added. Fig. 3 shows the comparison 
between an unfiltered record (3. a), with its corresponding group velocity (3.c), 
and filtered record (3.b) of a teleseismic event. 

The time-frequency filter efficiently reduces the influence of noise, body waves 
and higher mode Rayleigh waves. It also makes it possible to identify and 
exclude events whose fundamental mode of Rayleigh waves is not well sepa- 
rated from other waves. 10 events were rejected during this stage. A further 
12 events were excluded during the array analysis, leaving 26 events. Table 1 
contains the final event list, and figure 2 shows their distribution. The number 
of teleseismic events was too small to ensure a correct back-azimuth distri- 
bution, but there were events from all quadrants. The regional events were 
mainly located in the Pacific Coast subduction and the Caribbean Islands, 
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ensuring a back-azimuth range between N126° (South East) to N273° (West). 



3 Methodology 



To measure phase velocities, we follow iPedersen et al.l (120031 ) . In this method 



the phase velocity at a given frequency is obtained in two steps. Firstly, each 
event is analysed independently. In this step, the phase of the Wiener filtering 
W(f) is transformed into time delays At between each pair of stations using 
At = 0/(271"/). The Wiener filtering in the frequency domain that we use is 
given by: 

^ Sxy * Han(f)e^ 



S X x * Han(f)^S Y Y * Han(f) 

Sxx, Syy and Sxy are respectively the Fourier transforms of the autocorrela- 
tions and intercorrelation of the two signals , Han(f) is the Fourier transform 
of the Hanning fonction han(t) and t is the delay of the intercorrelation 
peak. * represents the convolution product. A grid-search on velocity and 
back-azimuth is applied to find the best fitting plane wave that would ex- 
plain the observed time delays. The fit is calculated with the LI norm, i.e the 
average absolute difference between observed and predicted time delays. The 
knowledge of the back-azimuth makes it possible to subsequently calculate the 
distance between each pair of stations projected onto the slowness vector. In 
this way each event yields a series of (d istance, delay) point s. 



Secon dly, a bootstrap process is applied (ISchorlemmer et al. 



2003 



Efron and Tibshiranil . 



19961 ): 500 bootstrap samples are created by resampling the 26 events of the 
data set. For each bootstrap sample, the phase velocity is calculated as the 
inverse of the slope of the best fitting line through all the (distance, delay) 
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points. The LI norm is also used here to estimate the fit between observa- 
tions and predictions. The points are associated with weighting which reflects 
how well the data fitted the assumption of a plane wave in the first step of the 
analysis. The final phase velocity and the associated uncertainity are obtained 
as the average and standard deviation over the 500 samples. 
The advantages of this method are 1) stabilization of delay measurements 
through Wiener filtering; 2) stabilization of back-azimuths and phase veloci- 
ties through the use of the LI norm; 3) weighting of the events in the final 
phase velocity calculation according to the quality of the back-azimuth esti- 
mate; 4) estimation of r ealistic error bars on the final dispersion curve. For 



more details, we refer to iPedersen et al.l (120031 ) 



The last part of this analysis consists in inverti ng the dispers i on cu rves. We 



used the two-step inversion methods proposed by 



Shapiro et al. 



(119971 ). Firstly, 



the average d ispersion curve w as inverted using a linearized, classical disper- 



sion scheme (iHerrmann 



19871 ) to find a simple shear wave model which fit- 
ted the dispersion curve. We then used this model for a stochastic nonlinear 
Monte-Carlo inversion. Interface depths and shear-wave velocities were ran- 
domly changed into a new model which was kept and used in the next iteration 
if it fitted the dispersion curve within the error bars. This second step was 
repeated 5000 times. We eliminated unrealistic models by applying loose con- 
straints on Moho depth (between 40 and 50 km depth) and on the S-wave 
velocity near the surface in agreement to the existing models (between 1.5 
and 3 km/s). We finally calculated the average model and we verified that the 
corresponding dispersion curve fitted within the error bars of the observed one. 
The error bars of the final models were computed as the standard deviation 
of all the acceptable models. Quality factors and P-wave velocities were kept 
constant during the inversion as their influence was significantly smaller than 
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the error bars of the dispersion curve. 



4 Results 



To detect differences between the crust under the Popocatepetl and the stan- 
dard crust of the MVB, one can compare the equivalent shear wave velocities 
profiles. As surface wave inversions are non-unique it is however useful to also 
compare the dispersion curves which correspond t o the existing mode 



Cruz-Atienza et al. 



s. 



(j200lh 



The models that we compare with were from: 1 
who obtained their model through inversion of receiver functions using four 
teleseismic events from South America a t station PPIG (locat ed 5 km north 
of the Popocatepetl crater, see fig.l); 2) 



Campillo et al. 



(119961 ) who inversed 



the gro up velocities o f loca l events between the Guerrero Coast and Mexico 



City; 3) 



Valdes et al.l (J1986J) whose model is the result of a seismic refraction 



study i n Oaxaca. We recalcu lated the phase velocities corresponding to these 



models. 



Shapiro et al. 



(119971 ) detected lateral variations of uppermost crustal 
structure within the MVB using surface wave group velocities. Due to the lim- 
ited depth penetration in their study (10 km), their models are not included 



in our figures, but will be integrated in the discussion of the results. 



4-1 Full array 

We firstly used all the stations and the 26 events to measure the 'overall disper- 
sion curve', i.e. the average dispersion curve within the full array. The phase 
velocities were unstable above 35 s period and were not used in the inversions. 
In Figure 4. a we compare our dispersion curve with the ones corresponding to 
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the earth model derived by 



and 



Valdes et al. 



Campillo et al. 



(1996) 



Cruz-Atienza et al. 



19861 ). For periods longer than 8 s, the 



Campillo et al 



(2001) 



(1996) 



curve is similar to ours. For short periods, the ve 



period, similarly to the 



Cruz-Atienza et al. 



ocities increase rapidly with 



(l200lh curve. 



We verified that our inversions of the observed phase velocities were indepen- 
dent of which of the three reference models (see fig 4.b) was used as starting 



mode 



Campillo et al. 



. The results shown here are obtained by using the one of 
( 119961 ) as it has the advantage of fitting our data well and it only has four 
layers. The latter is important to allow for an efficient exploration of the pa- 
rameter space in the Monte Carlo inversion. 

Our preferred model (fig. 4.b) shows low shear velocities (2.2 km/s) between 
the surface and 3 km depth, overlying a layer with velocities increasing slowly 
from 3.4 to 3.7 km/s between 6 to 20 km depth. The transition between the 
two layers may be either a strong gradient or a sharp interface. The lower 
part of the crust has a constant velocity of 3.75 km/s down to Moho which 
is located at 45 km depth. The velocity below Moho is approximatively 4.3 
km/s. The lower crust and upper mantle velocity as well as the Moho depth 
are not well resolved due to trade-off between these parameters and because 
of a maximal period of 35s. 

The boundary depths that we obtained are howev er consistent with existing 



models, in particular with 



Campillo et al. 



(119961 ). Our near surface veloci- 



ties are however sig nificantly lower a nd ou r upper crustal velocities slightly 



higher than those of 



Campillo et al 



(119961 ). while the two models are virtu- 



ally almost identical in the lower crust. The low velocity of the surface layer is 
relatively well constrained, however we can not exclude that the layer would 
be slightly thinner with s lightly lower velocity. This layer, also identified by 



Cruz-Atienza et al. 



(]200ll ). can be associated with the poorly consolidated ma- 
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terials of the volcano cone. 



Shapiro et al. 



(119971 ) find that the velocities in the 



upper 2 km are low beneath the southern part of the MVB where the volcanic 
activity is recent as compared to the northern part. Our results imply that 
the overall crustal structure below Popocatepetl is not significantly different 
from that of the MVB. 

We verified whether an 6-layers ini tial model with a Low ve locity zone be- 



tween 6 and 10 km inspired by the 



Cruz-Atienza et al. 



(1200 lh model, would 



yield a significantly different result. The resulting model is not different from 
our prefered model (fig 4.b), in particular there is no significant Low Velocity 
Zone. We do not see any indication that the Low Velocity Zone observed by 



Cruz-Atienza et al. 



(1200 ll ) is a general feature of the volcano. 



4-2 Sub-arrays 



To investigate lateral variations within the area, we divided the array into 
sub-arrays for which we calculated dispersion curves independently. The use 
of sub-arrays was particularly difficult as these arrays were composed of only 
three stations, so technical problems at any of the relevant stations would 
render the analysis impossible. It was possible to measure dispersion curves 
for the Southern (South sub-array: FPC, FPP, FPX) and the western sub- 
array (West sub-array: FPA, FPP, FPX). The dispersion curves for these two 
sub- arrays are shown in figure 5. At the largest period the analysis is mainly 
based on teleseismic events, out of which only 3 or 4 events were avalaible for 
the sub-array analysis. The phase velocity error bars are consequently large 
at long period (> 25 s for the South sub-array and > 15 s for the West sub- 
array) . 
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The dispersion curves for the South sub-array ( located around the active 
crater) and the West sub-array was respectively measured with 13 and 14 
events (see table 1). For the West sub-array, individual dispersion curves show 
strong oscillations between 6 and 12 s, particularly for events coming from the 
South or the East. These oscillations, probably due to local diffraction result 
in large error bar for the final dispersion curve. However, the dispersion curves 
beneath the two sub-arrays are not significantly different from the overall dis- 
persion curve. 

The phase velocities obtained with events for which the surface waves propa- 
gated through the volcano before encountering the array are more fluctuating 
than those obtained with other events. As the majority of the events are lo- 
cated South and South- West of the array, the individual dispersion curves are 
more fluctuant with the period at the North and East sub-arrays (composed 
of stations FPA, FPP, FMI and FPC) than for the other sub-arrays. It was 
consequently not possible to calculate a stable dispersion curve for these two 
arrays. 

The starting model for the inversion for the sub-arrays South and West is the 
model found with the full array, approximated by four layers. The estimated 
velocities are not significantly different from those of the overall model. Nev- 
ertheless, for the South array velocities are slightly smaller between 6 and 10 
km depth, and the surface velocities are higher. The differences are however 
smaller than the error bars of the overall model. 
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5 Discussion and Conclusion 



The average crust beneath the Popocatepetl volcano appears to be similar to 
the MVB crust. There is therefore no indication of large scale crustal anoma- 
lies associated with Popocatepetl as compared to the MVB. We do however 
confirm that the lower crust in the area is likely to be associated with rela- 
tively low shear wave velocities (3.75 km/s). Close to the surface, the velocity 
is approximatively 2.2 km/s over at least a depth of 3 km. It probably corre- 
sponds to the poorly consolidated material of the cone (such as volcanic slags 
and ash and py roclastic depos i ts) oy erlying the 2 km-thick volcanic layer of 



the MVB crust (jShapiro et al. 



1993). 



We speculate that the oscillations observed for the sub-arrays between 6 and 
12 s periods is associated with diffraction by lateral heterogeneity at 5-10 km 
depth as this period range corresponds to wavelengths between 16 and 36 km. 
To obtain strong diffraction, the heterogeneity must be of considerable size (i.e. 
of the order of the wavelength) , as surface waves are not strongly diffracted 



by many small heterogeneities (jChammas et al. 



20031 ). However, the lack of 



Low Velocity Zone turns down the hypothesis of a large continuous magma 
chamber. We speculate that either the interface located at 4 km depth in the 
average model fluctuate strongly, or that an abrupt lateral change takes place 
immediately beneath the central part of the volcano. The unresolved velocities 
at 5- 10s period at the West and South sub-arrays indicate that future arrays 
should be designed so as to give good constraints at 5-10 km depth. To obtain 
this, more stations and a larger recording period are necessary. 
More seismic events with a better azimutal distribution would improve the 
smoothing and the error bars of the dispersion curves and make it possible to 
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include receiver function analysis and coupled Rayleigh-Love inversion. 
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TABLE AND FIGURE CAPTIONS 



Table 1: Origine time, epicentre distance and back-azimuth of the events used 
for measuring the dispersion curves for the full array or with the sub-arrays 
(columns 5 and 6). 



Figure 1: (a) Location of the Popocatepetl volcano and (b) array geometry 



used i n the analysis. PPIG is a permanent station used by 



Cruz-Atienza et al 



( I200ll ) and is not used in this study. 



Figure 2: Location and azimuth distribution of (a) teleseismic and (b) local 
events used in the array analysis. 



Figure 3: Example of frequency-time filtering: a) Trace recorded at FPX, of 
the event at 03:37:42 GMT on november 03th 2002; b) Same trace after fil- 
tering; c) Group velocity of this event before filtering. 



Figure 4: a: Comparison of dispersion curves: 

1) Uncertainities of our observed dispersion curves (± 1 standard deviation, 
grey area) for the full array and 2) dispersion c urve calculat i ng wi th our av- 



era ge model (solid 



for 



Campillo et al. 



l ine) ; 



3) Dispersion curve for 



(119961); 5) Same for 



Valdes et a 



Cruz-Atienza et al. 



fj200lh . 



L986J); 4) Same 



Fig. 4. b: Comparison of crustal models: 

1) Error bars (± 1 standard deviation, grey ar ea) and 2) av e rage S-wave ve- 



locity model (solid line); 3) Crustal model for 



Valdes et al. 



(119861 ): 4) Same 
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for lCampillo et al.l (119961 ); 5) Same for ICruz-Atienza et al.l (120011 ). 



Figure 5: Dispersion curves of the full array (dotted line) and the two sub- 
arrays (solid line) with their uncertainities (grey area): a) South sub-array; b) 
West sub- array. In insert: sub-array geometry and back-azimuths of the events 
used. 
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